Part 2: Can more complex models wtih land cover and/or weather data improve on the best basic model?
Land Cover Model 1 - Urban Fabric Continuity Only
s1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = sun)
w1_lc <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + scale(uf_cont_destino) + scale(uf_cont_residencia), data = wed)
Land Cover Model 2 - All Land Cover Vars
lc2_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))
s2_lc <- lmer(formula = lc2_formula, data = sun)
w2_lc <- lmer(formula = lc2_formula, data = wed)
Weather Model
s3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = sun)
w3_wea <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip, data = wed)
Land Cover & Weather Model
all_formula <- as.formula(paste0("log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city) + temp + I(temp^2) + precip + ", paste0("scale(",names(sun)[12:23],")", collapse = " + ")))
s4_all <- lmer(formula = all_formula, data = sun)
w4_all <- lmer(formula = all_formula, data = wed)
Comparison of Complex Models to Basic Model
AIC(s_basic,s1_lc,s2_lc,s3_wea,s4_all)
## df AIC
## s_basic 7 552433.7
## s1_lc 9 552144.1
## s2_lc 19 537242.1
## s3_wea 10 551799.6
## s4_all 22 536568.5
AIC(w_basic,w1_lc,w2_lc,w3_wea,w4_all)
## df AIC
## w_basic 7 1144243
## w1_lc 9 1143452
## w2_lc 19 1129340
## w3_wea 10 1142632
## w4_all 22 1127414
Part 3: Leave Out City, Make Predictions, Compare to Observation
s_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(sun,city != "B"))
w_basic <- lmer(log(flujo) ~ scale(dens_destino) + scale(dens_residencia) + log(dist) + I(log(dist)^2) + (1 | city), data = filter(wed,city != "B"))
s_complex <- lmer(formula = all_formula, data = filter(sun,city != "B"))
w_complex <- lmer(formula = all_formula, data = filter(wed,city != "B"))
sun_pred <- sun %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(s_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(s_complex, ., allow.new.levels = T)))
wed_pred <- wed %>% filter(city == "B") %>% mutate(flujo_basic = exp(predict(w_basic, ., allow.new.levels = T)), flujo_complex = exp(predict(w_complex, ., allow.new.levels = T)))
map <- suppressMessages(st_read("~/GoogleDrive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp")) %>% filter(ID_GRUPO %in% unique(sun_pred$codigo_destino))
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
sun_res <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_des <- sun_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
sun_geo <- bind_rows(sun_res,sun_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
wed_res <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_residencia, lat = lat_residencia) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_des <- wed_pred %>% select(lon_residencia, lon_destino, lat_residencia, lat_destino, flujo, flujo_basic, flujo_complex) %>% bind_rows(setNames(.,c("lon_destino","lon_residencia","lat_destino","lat_residencia","flujo","flujo_basic","flujo_complex"))) %>% group_by(lon_residencia, lon_destino, lat_residencia, lat_destino) %>% summarise_all(sum) %>% ungroup() %>% mutate(id = row_number(), lon = lon_destino, lat = lat_destino) %>% select(id,lat,lon,flujo,flujo_basic,flujo_complex)
wed_geo <- bind_rows(wed_res,wed_des) %>% st_as_sf(coords = c("lon","lat"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84") %>% group_by(id) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% st_cast("LINESTRING") %>% mutate_if(is.numeric,log) %>% arrange(flujo)
p1 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = sun_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p2 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p3 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(sun_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Complex Model")
p4 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = wed_geo, aes(color = flujo)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Observed")
p5 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_basic), aes(color = flujo_basic)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Basic Gravity Model")
p6 <- ggplot() + geom_sf(data = map, color = "#404040", fill = "#101010") + geom_sf(data = arrange(wed_geo,flujo_complex), aes(color = flujo_complex)) + scale_colour_gradient2(low = "#191970", mid = "#00008B", high = "#ffff00", midpoint = 7.5) + theme_void() + theme(plot.background = element_rect(fill = "#101010"), legend.position = "none", plot.title = element_text(colour = "white", hjust = 0.5)) + labs(title = "Complex Model")
Barcelona
Sundays

Wednesdays

Madrid
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
Sundays

Wednesdays

Valencia
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
Sundays

Wednesdays

Sevilla
## Reading layer `celdas_marzo_2020' from data source
## `/Volumes/GoogleDrive/My Drive/Data/INE Mobility/Supporting Info/shapefiles_ reas_movilidad/Shapefile_celdas_marzo_2020/celdas_marzo_2020.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3214 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2021660 ymin: 3203377 xmax: 481778.5 ymax: 5433001
## Projected CRS: WGS 84 / Pseudo-Mercator
Sundays

Wednesdays
